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Abstract: We present a first-principles method to compute radiation 
properties of ultra-high quality factor photonic crystal cavities. Our 
Frequency-domain Approach for Radiation (FAR) can compute the far-field 
radiation pattern and quality factor of cavity modes ~ 100 times more 
rapidly than conventional finite-difference time domain calculations. It also 
provides a simple rule for engineering the cavity's far-field radiation pattern. 
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The high quality factors (Q) and small modal volumes of photonic crystal (PC) cavities 
make them ideally suited for applications requiring strong optical field enhancement, such as 
low-energy optical switching [ 1 ], strongly coupled cavity quantum electrodynamics (QED) (2) 
and harmonic generation [3 ]. Recent interest has focused on high-Q cavities where the far-field 
radiation pattern is engineered to emit vertically, enabling free-space mode excitation (U[5l[6l 
in cavity QED [ 5 ] and harmonic generation [7 ] experiments. 

Photonic crystal cavity design uses established theoretical ideas [9] [lOl [TT] Q2) to maxi- 
mize quality factors, in conjunction with finite difference time domain (FDTD) calculations to 
compute the cavity mode. Even with improvements in speed and accuracy |H3][T4|, time do- 
main calculations are by nature computationally intensive for the long life-times of ultra-high 
Q cavities, taking hours or days per design on a supercomputer. Optimizing both the quality 
factor and the radiation pattern can require the exploration of a large parameter space OO, 
further increasing the computation effort. These severe computational demands, and the inabil- 
ity of FDTD to provide insight into the underlying physics, point to the need for an alternative 
method. Here we provide such a method. Our first-principles Frequency-domain Approach for 
Radiation (FAR) is ~ 100 times more efficient than FDTD calculations since we do not compute 
radiative modes directly. It consists of two parts: we initially approximate the cavity mode as a 
bound mode after which the radiation is obtained using perturbation theory, thus avoiding the 
most time-consuming part of FDTD calculations. The FAR also provides a design strategy for 
achieving cavities with specified radiation patterns without requiring exhaustive simulations. 

We apply the FAR to double heterostructure photonic crystal cavities 1 10], which are formed 
by perturbing a photonic crystal waveguide (PCW) in a slab geometry. The two geometries are 
shown in Fig. [TJa)-(b); in the photosensitive cavity (Fig. [IJa)), the refractive index of a strip 
around the PCW (yellow shading) is uniformly increased by An p , as can be achieved in chalco- 
genide glass (ElIIS). In the fluid infiltrated cavity (Fig.[TJb)), the refractive index of the holes 
is increased by Ani in a strip-like region (red shading), typically by fluid infiltration |[T7l[T8ll . 
These two cavities are therefore complementary, i.e., in one only the background is perturbed, 
while in the other only the holes. We have found that considerable qualitative insight into the 
radiation pattern of the cavity mode can be obtained by examining a single term in the equation 
that governs the radiation from the cavity. This term has the form A(r)D fl (r), where A(r) is 
associated with the perturbation that creates the cavity and D a (r) is the bound approximation 
for the cavity mode. This term is shown for a z = slice through the PC slab in Figs.[TJc)-(d). 
The perturbation term A(r) is only non-zero in the background for the photosensitive cavity 
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Fig. 1. Schematic of (a) photosensitive cavity with locally increased background index and 
of (b) fluid infiltrated cavity with increased hole index. A(r)D fl (r) (arb. units) for (c) pho- 
tosensitive cavity (d) fluid infiltrated cavity, both with a length of L = Ad. (e)-(f) Modulus 
of the Fourier transform of (c)-(d) respectively, with non-radiating components removed. 



(Fig.[TJc)) and only non-zero in the holes for the fluid infiltrated cavity (Fig.[IJd)). The Fourier 
components in the light cone of this product are peaked near the edge of the light cone for the 
photosensitive cavity (Fig.[TJe)) corresponding to radiation being directed towards the horizon. 
However, for the fluid infiltrated cavity, the Fourier transform is strongest near k x = k y = (Fig. 
[TJf)), and thus it predominantly radiates vertically. We later return to this insight and use it to 
provide a general design rule for engineering the radiation pattern of a cavity mode. 

Our theory uses a Hamiltonian formulation fT9ll to construct cavity modes by superposing a 
basis of bound PCW modes expressed in terms of the B(r) and D(r) fields, so any superposition 
is divergence-free. The Hamiltonian for a dielectric PC cavity with relative permittivity e(r) is 

Since we use PCW modes as a basis, it is convenient to define e(r) = e(r) + e(r), where e(r) 
is the permittivity of the PCW, and e(r) the small permittivity change that creates the cavity. 
We then expand the cavity mode using the normalized PCW modes fT9l below the light cone 



D(r,0= / dkJ^a k e- ia *B k (r)+c.c., (2) 

Abound V I 

where we only include modes of the even PCW band, i.e. those for which E y (r) is even in y. 
We can include more modes, but ultra-high Q cavity modes are typically gently confined and 
thus different bands couple weakly. 

Substituting ^ into ([T]) we obtain an approximation for the Hamiltonian of the PC cavity 



a\a k t (3) 

where y(r) = l/(2eb) [l/e(r) — l/e(r)], and we dropped non-rotating wave terms involving 
a\a^ k , and a k a k i. Diagonalizing M{ determines an eigenvalue, the energy HcOq of a photon in the 



M{ = J dkdk' hcQ k S(k-k')+hy/aj^ty J dry(r)Dl(r) T> k 



(r) 



cavity mode, while its eigenfunction vo(k) gives the cavity mode in the basis of PCW modes: 



D fl (r)=/ dkJ^vo(k)D k (t). 

Abound V I 



(4) 



We now have an approximate expression for the cavity mode in terms of a basis with Fourier 
components outside the light cone. The Fourier content within the light cone of the ultra-high Q 
factor cavities of interest here is small, and we have found that D a (r) is a good approximation 
for the shape of the cavity mode. Similarly, the eigenvalues of ^ approximate the real part of 
the frequency of the cavity mode well. We thus use D a (r) to find a first approximation for the 
polarization field P(r) within the light cone. 

The polarization field P(r) = £o [e(r) — l]E(r) of a mode with frequency CO satisfying the 
macroscopic Maxwell equations is also a solution to the integral equation 



P(r) = £q [e(r) - 1] J dr f G(r - r'; co)F(r f ) 



(5) 



where the Green tensor expresses the electric field at y' due to an oscillating polarization source 
at r. We use the formalism for layered media l20lL in which we deal with a sheet of polarization. 
Since we need to compute the out-of-plane (z-direction) radiation of a PC cavity, this formalism 
is particularly useful as it separates propagating modes in the z-direction, with |ic| 2 = t% + ky < 
from evanescent modes with | k\ 2 > k^, where ko = COo/c. 

Defining T(r) = [e(r) — l]/e(r), the polarisation field of the cavity is approximated by 
P a (r) = r(r)D <2 (r) = (r(r) +f (r))D <2 (r), where again the over-bar denotes a quantity for the 
PCW and the tilde denotes the perturbation creating the cavity. Since D a (r) has no Fourier 
components within the light cone, neither does r(r)D fl (r), T(r) being periodic with the lat- 
tice. However, f (r)D <2 (r) does have components within the light cone, providing a starting 
point for calculating the radiative polarization. We relate the actual polarization field of the 
cavity mode P(r) to P fl (r) by writing P(r) = F a (r) +P c (r), whose radiative components are 
p rad (r) = f (r)D a (r) -f P£ ad (r), where P c (r) is the correction to the polarization field, while 
rad refers only to Fourier components in the light cone. We write the complex cavity mode 
frequency CO as CO = (Oq + CO. We next perform a Taylor expansion about COo of the Green func- 
tion, substitute into and use the fact that J*™ d (r) and the variables with tildes are small. 
After some manipulation and keeping only terms with Fourier components in the light cone, 
we obtain a first order expression for P rad 



PfV) -S)(e(r) - 1) y > Jr / G(r-r / ;o) )Pi ad (r / ) = A(r)D» 



f(r) g « f « 



e(r) 



D"(r), 



(6) 

where the driving term, which contains information about the cavity via the parameters with a 
tilde, couples to Fourier components inside the light cone. As discussed earlier in this paper, we 
have found that A(r)D fl (r) gives good qualitative insight into the far-field radiation. In general 
though, Eq. ^ is a Fredholm integral equation of the second kind, in which the Green function 
ensures a self-consistent interaction between the dipoles. 

By solving ([6]), we obtain the full quantitative radiative polarization components of the cavity 
mode, from which the far-field radiation can be determined using the Green function in (J5J. In 
the far-field we write the electric field as Ef ar (r) = e±(k)e lk ° r /r, where k = kor- (xx + yy), 
with, above (+) and below (— ) the slab, 



4(K) = i s • J dzdRe- iK - R e Tiwz F\ ad (R,z) (7) 
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Fig. 2. Quality factor versus cavity length for (a) the photosensitive cavity (Fig.[TJa)); (b), 
the fluid infiltrated cavity (Fig.[TJb)). Red symbols are computed using the FAR method, 
while blue ones are computed by FDTD. 



for s polarization, and with a similar expression for p polarization. Here R = (x,y) and w = 
I*! 2 ) 1 / 2 . Equation {7} is thus a planar (x and y) Fourier Transform, integrated over the 
thickness of the slab (z) with appropriate phases. Each (k Xj k y ) of the polarization field inside 
the light cone corresponds to a unique far-field direction. The far-field electric field gives the 
Poynting vector, and therefore the quality factor of the cavity mode can be computed. 

To obtain numerical solutions to ([6]), we further assume that Fourier components inside the 
light cone do not couple to those outside the light cone. Since inside the light cone k XJ are small, 
this lets us use a coarse discretization in x and y, reducing the size of the problem. By using 
an efficient iterative bi-conjugate gradient method |ED[22l, Eq. ([6]) can be solved to within a 
tolerance of 10 -5 in 20 — 100 iterations, each of which take less than 10 seconds. Our MATLAB 
code typically solves ^ and ^ in under 15 minutes on a work station. In contrast, the FDTD 
calculations for each point in Fig.[2]took tens of hours on a 32 core cluster. 

In our simulations for the photosensitive cavity (Fig.[TJa)), we take aWl PCW with back- 
ground index of = 2.7, slab thickness, t = O.ld and hole radius a = 0.3d, where d is the 
period and An p = 0.02,0.04. For the fluid infiltrated cavity (Fig.[T|b)) we use a W0. 98 silicon 
PCW (background index = 3.46), with slab thickness, t = 0.49 d, hole radius a = 0.26J, 
Ant = 0.2,0.4,0.6. In Fig. [2] we show the Q-factor versus cavity length calculated using the 
FAR method (red) and using FDTD (blue). In Fig.[2ja), which is for photosensitive cavities, the 
efficiency of our theory allows us to vary the cavity length continuously. This is impractical for 
FDTD calculations, so we only have results at even integer values of the cavity length and at 
some intervening points. The agreement between the results is excellent: the Q-f actors agree to 
within 30% (or their logarithms by 2%), making them suitable for examining trends in Q. The 
strong oscillations in Q correspond to a factor of 8. In Fig. |2jb), which is for fluid infiltrated 
cavities, we only calculated Q for even integer cavity lengths. The agreement for these cavities 
is good: the results have the same trends and never differ by more than a factor two. 

Having demonstrated the reliability of the FAR, we now exploit its semi-analytic nature to 
design desirable far-field radiation properties. Figure [3] shows good agreement between the far- 
field radiation patterns computed using the FAR (left) and FDTD (right), for photosensitive 
cavities (Figs.[3Ja),(b)) and fluid infiltrated cavities (Figs.|3jc),(d)) of different lengths L. Note 
that (i) the number of lobes in the radiation pattern increases as the cavity gets longer; and that 
(ii) as discussed earlier, photosensitive cavities radiate predominantly at large declination angles 
(6), while the fluid infiltrated cavities radiate mostly vertically. Both features can be explained 
by examining A(r)D a (r). The effect of A(r) on D a (r) is to introduce nodes and anti-nodes due 
to Fabry-Perot effects in the cavity. Point (ii) is more subtle: returning to Fig.[T] since the cavity 
modes are dielectric modes, for the photosensitive cavity the product A (r)D fl (r) (Fig.[T|c)) has 
the effect of merely introducing sidelobes in the Fourier transform of D <2 (r). The overlap of 



these sidelobes with the light cone (Fig.[T|e)), is dominated by (k x ,k y ) values at the edge of the 
light cone, maximizing the Q factor and leading to radiation at large declination angles. 




Fig. 3. Symmetric quadrants of far- field radiation (S r ) for (a),(b) cavities in Fig.[TJa) with 
An p = 0.02 and (c),(d) those in Fig.JTJb) with An/ = 0.2. Left frames are computed using 
the FAR method while right frames are computed using FDTD. Colors are as in Fig.fTJd). 
Angles and 6 are azimuthal and declination angles respectively. 



For fluid infiltrated cavities, A(r)D <2 (r) (Fig. [TJd)), is nonzero only inside holes. Its Fourier 
transform within the light cone (Fig.[TJf)) peaks at the origin, because the cavity length is such 
that A(r)D fl (r) has a strong non-zero DC Fourier component. This is clear from the fields in 
the holes in Fig.[TJd): four holes have strong positive fields and only two have strong negative 
fields because the cavity mode is dominated by the Bloch mode at kd = n, which changes 
sign each period. This cavity therefore radiates mostly vertically. Since examining A(r)D <2 (r) 
is sufficient for qualitative insight into the far-field, the requirement for steering the radiation 
of cavity modes is thus simple: construct a perturbation such that A(r)D fl (r) has a Fourier 
transform which peaks at (k x ,k y ) values corresponding to the desired direction. 

Similar arguments can be used to explain the variations in Q observed in Fig.[2ja): the cavity 
mode is a superposition of Bloch functions centred about the Brillouin zone edge (kd = it). It 
is thus not surprising that the period of the oscillations in Fig.[2ja) corresponds to the period of 
the central Bloch function. The details of the oscillations in Q depend on the superposition of 
the Bloch modes in A(r)D a (r) overlapping with the light cone. 

We have presented a frequency-domain approach for radiation (FAR) that allows the effi- 
cient calculation of the radiative properties of ultra-high Q PC cavities. Both Q-f actors and the 
radiation patterns are in good agreement with fully numerical FDTD calculations. The orders- 
of-magnitude improvement in computation speed will enable the application of powerful opti- 
mization algorithms, potentially transforming PC cavity design. The FAR lets us directly predict 
the radiation pattern through its link to the cavity's refractive index. Although we applied the 
theory to cavities created by refractive index changes, extensions allow the treatment of other 
cavity types, created, for example, by shifting inclusions [23 ] or stretching the lattice ifTOl . 
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